Quantitative Susceptibility Map Reconstruction with Magnitude Prior

INTRODUCTION: Quantitative Susceptibility Mapping (QSM) aims to quantify tissue magnetic susceptibility with applications such as tissue contrast enhancement [1], venous blood oxygenation [2], and iron quantification [3]. The magnetic susceptibility χ maps to the observed phase shift via a well-understood transformation, but the inverse problem, i.e. estimation of χ from phase, is ill posed due to zeros on a conical surface in the Fourier space of the forward transform; hence χ inversion benefits from additional regularization [4]. Here we propose enhanced regularization for χ inversion by incorporation of magnitude priors. Since the data acquisition step for QSM yields both phase and magnitude data, the inverse problem can be better conditioned if the magnitude is incorporated as a prior. By encoding spatial priors derived from a magnitude image into an l1 regularization scheme via the Focal Underdetermined System Solver (FOCUSS) algorithm [5], we report high quality QSM on a numerical phantom with four-fold improvement in RMSE when magnitude priors were applied. We also demonstrate the application of the method on in-vivo data at 7T. THEORY: The system of linear equations δ=FDFχ defines our ill-posed deconvolution task where δ=φ/(γ·TE·B0) is the normalized field map, D is the susceptibility kernel in k-space, F is the Fourier transform operator, χ is the susceptibility vector, φ is the unwrapped phase, γ is the gyromagnetic ratio, TE is the echo time and B0 is the main field strength. We assume that χ shares tissue contrast boundaries with the magnitude image, and is therefore expected to have similarly sparse spatial gradients as the magnitude; this prior knowledge can be imposed on the reconstruction via the regularized FOCUSS algorithm [6]. Letting ∂x χ denote the spatial gradient along x, the k step of the iterative algorithm is as follows: set Wk+1=diag(|∂x χk|), solve qk+1=argminq ||VxFδ−DFWpriorWk+1q||2 + λ||q||2, and update ∂xχk+1=WpriorWk+1qk+1. Here, Vx is a diagonal matrix that acts as gradient operator in k-space due to Vx(ω,ω)=(1−e) where n is the matrix size along x. The diagonal weighting matrix Wprior=diag(|∂xm|) is generated from the magnitude image m to express our prior belief that the magnitude and susceptibility images share similar gradients. This is seen when the Least Squares (LS) solution step is expressed as a function of ∂x χ as follows: ∂xχk+1=argmin∂χ ||VxFδ−DF∂χ||2 + λ||WpriorWk+1∂χ||2. Here, when the magnitude gradient ∂xmi at voxel i is small, Wprior(i,i) will be large and penalize ∂xχi more. After obtaining the susceptibility gradients in three dimensions, we estimate χ by solving a LS problem: χ=argminθ Σr=x,y,z ||∂rχ−∂rθ||2 + β||δ−FDFθ||2 where we used β=1 in our experiments. METHODS: The numerical susceptibility phantom (x×y×z=128×128×32) in Fig. 1 contains three compartments: a rectangular prism (χ=1ppm), a cylinder (χ=0.047ppm) simulating graywhite matter susceptibility difference [1], and a 2-pixel wide vessel (χ=0.4ppm). The vessel has three segments; i) along B0 (z-direction), ii) in-plane part, and iii) a 35° slanted segment, which is perpendicular to the magic angle of 55° and therefore poses the most challenging inversion geometry. We also created a magnitude image with shared boundaries, but with different compartment intensities. Starting from the true susceptibility, we forward simulated the field map by convolution and corrupted it with complex valued Gaussian noise so that the noisy field map had 17.9% normalized root mean squared error (RMSE) relative to the noise-free case. We tested the FOCUSS algorithm without a prior (by setting Wprior = I) and with magnitude prior and used an optimal λ setting (λ=10 without prior and λ=10 with prior) to reconstruct the susceptibility from the noisy field map. Second, we tested the FOCUSS on in-vivo data. At 7T, a 3D GRE sequence was used to acquire axial images with 0.33 mm in-plane resolution, 1.0 mm slice thickness and FOV of 192×168×64 mm for a TE of 10 ms on a young (26 years, female), healthy subject. After high-pass filtering the phase with a Hanning filter of size 64×64, the susceptibility distribution was reconstructed from the field map using the FOCUSS algorithm (λ=10) by using the magnitude as a prior. We report susceptibility differences Δχ = χvessel – χtissue for a selected vessel and Δχ = χputamen – χtissue for the putamen by manually generating interior and surrounding tissue masks for averaging. RESULTS: For the numerical phantom, susceptibility reconstruction without a prior resulted in 5.2% RMSE whereas using the magnitude prior yielded 1.3% error (Fig. 1). We note that the slanted segment of the vessel is almost invisible in the field map due to the ill-posed kernel, but FOCUSS with magnitude prior successfully recovered this segment. Fig. 2 depicts 7T QSM results obtained after taking maximum intensity projection (MIP) over 4 slices that contain the vessel of interest. We computed χvessel by taking the mean of the MIP image inside the vessel ROI (Fig. 2b) and χtissue by taking the mean susceptibility of pixels inside the tissue ROI across 4 slices. In this case, Δχ was estimated to be 0.34ppm. Fig. 3 presents QSM results for the putamen obtained by taking the average susceptibility over 6 slices. After computing average χputamen and χtissue inside the ROIs in Fig. 3b-c, Δχ was estimated to be 0.011ppm. Our result falls within the range of putamen susceptibility values (0 to 0.054ppm) reported in [3] for subjects between 20 and 30 years. CONCLUSION: By making use of magnitude information to add spatial priors to l1 regularization, we demonstrate high quality QSM on numerical and in-vivo data. In addition to estimating venous oxygenation, the algorithm can be used for quantification of susceptibility inside iron-rich brain structures. REFERENCES:[1] Duyn JH et al., PNAS 2007;104(28):1179611801 [2] Fan AP et al., ISMRM 2010;693 [3] Liu T et al., ISMRM 2010;4364 [4] Liu J et al., ISMRM 2010;4996 [5] Gorodnitsky IF et al., IEEE T Signal Proces 1997;45(3):600-616 [6] Cotter SF et al., IEEE T Signal Proces 2005;53(7):2477-2488 −1 −1